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Abstract. In this paper a macroscopic model of tumor cord growth is developed, relying on 
the mathematical theory of deformable porous media. Tumor is modeled as a saturated mixture 
of proliferating cells, extracellular fluid and extracellular matrix, that occupies a spatial region 
close to a blood vessel whence cells get the nutrient needed for their vital functions. Growth 
of tumor cells takes place within a healthy host tissue, which is in turn modeled as a saturated 
mixture of non-proliferating cells. Interactions between these two regions are accounted for as an 
essential mechanism for the growth of the tumor mass. By weakening the role of the extracellular 
matrix, which is regarded as a rigid non-remodeling scaffold, a system of two partial differential 
equations is derived, describing the evolution of the cell volume ratio coupled to the dynamics 
of the nutrient, whose higher and lower concentration levels determine proliferation or death 
of tumor cells, respectively. Numerical simulations of a reference two-dimensional problem are 
shown and commented, and a qualitative mathematical analysis of some of its key issues is 
proposed. 



1. Introduction 

Mathematical modeling of tumor growth dates back to many years, as the considerable liter- 
ature on the subject demonstrates. Recent surveys by Araujo and McElwain [3] and by Preziosi 
and Graziano |19j highlight the increasing interest by biologists and mathematicians toward the 
qualitative and quantitative comprehension of this complex phenomenon, with the common aim 
of devising suitable descriptive tools to predict the behavior of the system and to virtually test 
the effect of different, possibly new therapies. 

As a matter of fact, it is all but an easy task to conceive a mathematical model capable 
to take into account the wide variety of biomechanical and biochemical factors, in most part 
intimately correlated to each other, that contribute to the overall dynamics of the system. This 
is particularly evident if one considers that tumor growth is essentially a multiscale phenomenon, 
in which subcellular, cellular, and tissue levels are strongly interconnected in terms of cause and 
effect. In many cases, the reason for some macroscopic outcome relies on the activation of specific 
internal mechanisms to the cells. These usually trigger particular collective behaviors of groups of 
cells, resulting finally in an observable effect at the macroscopic scale. In addition, when dealing 
with living matter, as cells and tissues of human body are, classical laws of Newtonian mechanics 
are not in principle the only rules that the evolution of the system obeys to, due to some sort of 
intelligence that makes cells able to self-organize according to possible changes in their state or in 
the environment they live in. 

Without claiming to be thorough, and referring instead to the above-cited papers for a compre- 
hensive review on the state of the art, we briefly outline here some basic facts about mathematical 
modeling of tumor growth at the macroscopic scale, in order to sketch the context which the 
present work flts in. 

A first thread of models of tumor growth to be mentioned (see e.g., Byrne [12] and the main 
references listed therein) is based upon a particular geometry of the tumor mass, the so-called 
spheroid, which essentially consists in a spherically-shaped three-dimensional aggregate of cells, 
primarily produced in vitro for experimental purposes. In this specific configuration, growth is 
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addressed by focusing on the time evolution of the external radius R(t) of the spheroid. An integro- 
difFerential equation is obtained by equating the time variation of the volume of the spheroid to 
the overall cell proliferation, under the assumption that the latter is somehow related to the 
concentration of a chemical, usually oxygen, which provides the cells with the nutrient needed for 
their vital functions. This way it is unnecessary to explicitly track the evolution of the cell density, 
and the model is closed by simply adding a stand-alone reaction-diffusion equation for the oxygen 
concentration, however posed in a domain which changes in time according to the evolution of 
R{t). The resulting mathematical problem is an integro-differential free boundary problem, which 
has been proved to admit solutions with suitable regularity properties and to predict, in some 
proper ranges of the parameters, the evolution of the system toward a steady state characterized 
by a finite steady growth size of the spheroid (Bueno et al. [T^, Friedman and Reitich [16]). 

A slight variant of this framework entails the introduction of two regions in the spheroid, namely 
an outer viable zone in which cells proliferate and a central necrotic core in which cells starve and 
die due to an insufficient delivery of oxygen from the periphery of the tumor. In this case, specific 
reaction-diffusion equations for the nutrient concentration are set up in either region, accounting 
for the fact that living cells take up oxygen from the environment, while necrotic cells only receive 
a small amount of oxygen by means of diffusion. In addition, an integro-differential equation is 
written for the evolution of the outer proliferating shell, in the same spirit as the corresponding 
equation of the previous model. The resulting problem describes once again the evolution of the 
nutrient concentration across the spheroid, but it is characterized by two free boundaries delimiting 
the necrotic core and the external periphery of the spheroid, respectively. From the mathematical 
point of view, the problem has been extensively studied by Cui and Friedman [T3], who addressed 
the existence of a stationary solution and the convergence to it of the time dependent model under 
suitable assumptions on the parameters. 

An alternative class of mathematical models of tumor growth relies instead on the theory of 
deformable porous media. In such a context, the tumor is regarded as a mixture of various mutually 
interacting components, which obey standard coupled mass and momentum balances. One of the 
main novelties with respect to the above-described spheroid models is that now no particular 
geometry is in principle required to write the equations, so that the global conformation of the 
tumor is in turn an unknown of the problem which can be duly studied by the model itself. 
However, the assumption of specific geometrical properties, like e.g., spherical or axial symmetry, 
or even the reference to one-dimensional settings are customary in the literature also in this case, 
in order to address simplified explorative models. 

Mixture-based models usually include one or more state variables specifically devoted to track 
the evolution in time and space of the cell population. In the early models of spheroids this 
was not a major issue since, as recalled above, the growth of the tumor was described at an 
essentially geometrical level by invoking a volume balance. Ambrosi and Preziosi [2 pointed out 
that this amounts to assuming a constant cell volume ratio within the tumor, so that the growth 
of the spheroid would be essentially dictated by the necessity to preserve such a constraint in 
spite of cell proliferation. Therefore, those models can be viewed as particular approximations of 
a more general multiphase setting. In many cases, the tumor is modeled as a biphasic mixture 
consisting of cells and extracellular material. Frequently, the former are further divided into 
several subpopulations of proliferating, quiescent, and necrotic cells, with possible interchanges 
among the populations determined by the availability of some nutrients that cells need in order to 
carry out their vital functions. The extracellular material is commonly identified with extracellular 
fiuid or extracellular matrix (stroma). It fills the interstices among the cells, so that no voids are 
left within the mixture and the latter can be regarded as saturated. In some models of vascular 
tumors, the extracellular material also includes a distribution of blood vessels mimicking the co- 
opted vasculature under angiogenic effects (Breward et al. [5]). 

In the theory of deformable porous media, the velocities of the constituents of the mixture are 
determined from suitable momentum equations accounting for the internal stress of each phase, 
as well as for the mutual external interactions among different phases. Therefore, as usual in 
continuum mechanics, suitable constitutive relations for the components of the tumor have to be 
specified, in order to render at the tissue level the proper mechanical behavior. It is plain that 
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this poses the major modehng problems, because hving tissues can be assimilated to classical 
materials only on first approximation. For a detailed study of the mechanics involved in tumor 
growth we refer the interested reader to Ambrosi and MoUica jlj. Here we simply remark that 
in some particular geometries mass balance equations may be sufficient to determine, at a purely 
kinematic level, the velocities of the components, at least under suitably reinforced assumptions 
on the saturation of the mixture (Bertuzzi et al. [5]). 

A rigorous derivation of model equations for avascular tumor growth according to the principles 
of the theory of mixtures can be found in Byrne and Preziosi and in Breward et al. [5]. 
Following analogous theoretical guidelines, in this paper we are concerned instead with vascular 
tumor structures, that develop mainly in the axial direction along the wall of a blood vessel. Due 
to the particular shape they take, they are named tumor cords. For such tumors no particular form 
of angiogenesis is needed, as they directly catch the nutrient by diffusion from the vessel which 
they grow around. Mathematical modeling of tumor cord growth has already been addressed by 
a number of papers by Bertuzzi and coworkers [SlEJ[7j. In particular, they use a mixture theory 
approach supplemented by proper kinematic constraints that allow to express the velocity of the 
components of the mixture without invoking any momentum balance. In addition, they assume 
cylindrical symmetry of the cords around the blood vessels, further neglecting axial variations of 
the state variables or performing suitable averages in the axial direction when needed, so as to 
reduce the domain of the problem to an annulus representing the two-dimensional cross-section of 
the cord. They also add a sophisticated kinematic constraint to describe the evolution of an inner 
viable core in interaction with an outer necrotic ring, which are modeled as two distinct zones 
by means of specific equations for each of them. Finally, they propose a detailed mathematical 
analysis of both the stationary and the evolution problems generated by their model. Conversely, 
we are interested in modeling tumor cord growth starting from sufficiently general principles of the 
theory of mixtures, so as to provide a model able in principle to face many different physical and 
geometrical configurations. A basic (i.e., minimal) version of the model is developed and studied 
here. For further extensions and improvements, as well as for more applications to different 
systems, we refer to Preziosi and Tosin [20]. 

The paper is organized into six more sections that follow this Introduction. In Sect. [2] the 
equations of the model and the proper boundary conditions are derived, then in Sect, p] the 
mathematical problem is stated in non-dimensional form. In Sect. [4] numerical simulations are 
shown and commented, with the aim of testing the ability of the model to reproduce some basic 
relevant features of the system. Sections |5] |6] address the qualitative mathematical analysis of 
some key issues of the problem, and Sect. [7] finally draws conclusions and briefly sketches research 
perspectives. 



2. The mathematical model 

In this section we propose a mathematical model for the description of the growth of a tumor 
mass along a source of nutrient, working in the framework of the theory of mixtures. In doing this, 
one of our main goals will be the minimality of the model, meaning that we will constantly aim 
at taking into account only those mechano-chemical mechanisms strictly relevant for an essential 
representation of the macroscopic phenomenon. Of course, we are aware that in this way many 
interesting aspects, whose inclusion would make the model closer to physical reality and perhaps 
less academic, are left out of the study. However, we believe that a minimal model along with 
a related qualitative mathematical analysis are necessary starting points before tackling more 
complicated situations, and that such a minimal model may be eventually used as the core of 
future more accurate models. 

For the sake of definiteness, we will deduce the model having in mind the two-dimensional 
geometry depicted in Fig. [l] Actually, since the model is intrinsically multidimensional, any 
further extension to different geometries and possibly to higher dimensions is mostly technical, 
requiring in principle the same ideas up to a suitable reinterpretation of the symbolic form of the 
differential equations that we will write. 
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Figure 1. The two-dimensional domain of the problem. A spatial region Q is 
divided into two subregions il and fi'^, representing the tumor cord and the host 
tissue respectively, separated by a free boundary S. The blood vessel coincides 
with the horizontal z axis. 

Let us then consider a spatial region Q C occupied by a mixture of cells, extracellular fluid, 
and extracellular matrix (ECM), that we suppose saturated. Denoting by 0c, '^m the respective 
volume ratios of the constituents, this condition is expressed by 

(1) 0c + (/if + 0m = 1 in Q. 

The region Q is internally divided into two subregions fi, fi'^ separated by a one-dimensional 
manifold 5, so that Q = 51 U fi"^ and n fi"^ = S*. We assume that represents the growing 
tumor cord and fi"^ the surrounding host tissue, while S plays in this context the role of a free 
boundary delimiting the cord. As we will see, the model guarantees that tumor cells located in Q 
and normal cells located in 17"^ never mix, hence the sole volume ratio 0c is sufficient to track the 
evolution of the whole cell population. The blood vessel around which the tumor cord develops is 
schematized by the horizontal z axis, and its presence will be incorporated into the equations of 
the model in terms of boundary conditions. This is possible, and is indeed useful, because we are 
not interested in the interactions of the vessel with the tumor, so that its geometry and mechanics 
can be duly ignored. 

2.1. Momentum equations. It is assumed that normal cells and tumor cells differ only in that 
the latter undergo an unregulated proliferation, while mechanical properties are the same for both. 
Hence, one can write a unique momentum equation of the form 



• p is introduced as a Lagrange multiplier to satisfy the saturation constraint (jlj, and is 
then interpreted as the interstitial pressure of the extracellular fluid; 

• Tc is the excess stress tensor of the cellular matter, accounting for cell-cell internal stress; 

• mc is the resultant of the actions on the cellular matter due to the interactions with the 
extracellular fluid and the extracellular matrix. 

Notice that in Eq. ([2| inertia is neglected, in view of the relatively slow dynamics of the cells 
during their growth process. 

Assuming that cells behave like an elastic fluid, we introduce a scalar, possibly nonlinear func- 
tion E = S(0c) such that 



(2) 



V • (0cTc) 4- 0c Vp = mc 



in Q 



where: 



(3) 



Tc - -S(0c)I, 



I being the identity matrix. Positive values of E for high 0c denote compression of the cells, while 
negative values for low 0c imply a packing tendency. In addition, we postulate the existence of a 
value 00 > such that S(0o) = 0, identifying a conflguration in which cells get in touch without 
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exerting any stress on each other. In the sequel, we will occasionally refer to 00 as the stress-free 
cell volume ratio (see Fig. |2|. 

Concerning the external actions, it is convenient to split nic into two contributions xiid, iHcm 
specifically related to the interactions of the cells with the fluid and the ECM. More in detail, we 
can imagine a simple viscous friction among the components of the mixture, that is 



where Vc, v^, denote the velocities of the constituents. In particular, in the classical Darcy-like 
theory of porous media it is customary to express the interaction with the fluid phase as 

(4) mrf= -^(vf-vc), 

with K > representing the permeability of the mixture. Regarding the interaction of the cells 
with the extracellular matrix, the theory is instead by far less classical and consolidated, and 
certainly viscous friction can be possibly accepted only as a flrst approximation. Indeed, cells 
are known to attach to the ECM via suitable adhesion sites, and to detach only in presence of 
sufficiently strong forces (Baumgartner et al. 4 , Canetta et al. [T3], Sun et al. [H]), therefore 
such an interaction is deflnitely not a pure sliding. However, to keep things simple we set 

(5) nicn = Ac„i(v„i - Vc) 

for a certain constant Ac„i > 0, referring the interested reader to Preziosi and Tosin [20J for a 
more detailed mathematical modeling of this term of the equations. 

Introducing Eqs. ^ into Eq. ^ we flnd, after some standard manipulations, 

(6) V(0cS(0c)) + 0cVp= ^(vf-Vc)+Ac„.(v,„-Vc) inQ. 

For the extracellular fluid, an equation similar to Eq. ([2| can be written: 

(7) - V • {^tTi) + 0fVp ^mi in Q 

with = m£c + nif„j, that is the sum of the interactions with the cells and the ECM, respectively. 
However, due to the interpretation of p provided above, the excess stress tensor is commonly 
neglected and the pressure is regarded as the main internal stress of the fluid, so that Eq. (|7| 
actually simplifies as 

(8) (j)i^p = nifc + in Q. 
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Relying again on a Darcy-like framework, we think of the interaction terms as proportional to 
the relative velocities of the interacting constituents, with specifically 



(9) 



due to an action-reaction principle, and 
(10) 



K 







to emphasize the fact that the interaction of the extracellular fluid with the extracellular matrix 
is actually negligible with respect to the corresponding one of the cells (cf. Eq. It might 

be inferred that a more accurate modeling of the latter may allow for a better expression of the 
former. However, such topics are beyond the scope of a minimal model of tumor growth and will 
therefore not be addressed in the present work. 

Putting Eqs. ([To| into Eq. (|8| yields the relation 



(11) 



which represents the classical form of the well-known Darcy's law. 

Finally, an equation similar to Eqs. ([2]), Q holds in principle also for the ECM: 



(12) 



in Q, 



with xiim = ni„jc + ^mi — —iHcm ~ m^m in view of the action-reaction principle. However, it 
requires to specify the excess stress tensor T^, which appears to be a very complicated object due 
to the fibrous nature of the extracellular matrix and to the continuous remodeling it undergoes 
when cells move within it. Therefore, wishing to focus mainly on the growth of the tumor cord in 
interaction with a source of nutrient, we choose to consider the ECM as a rigid non-remodeling 
scaffold in which cells live and grow and extracellular fluid flows. As a consequence, we replace 



Eq. ([12f by 

(13) v„j = in Q 

and we simply observe that plays now the role of a (tensor) Lagrangian multiplier to satisfy 



the constraint (13 1. In addition, we take the volume ratio 

(14) 0™ = l-0*, 
so that the saturation constraint ([T]) specializes as 

(15) 0c = 



as a constant, say 
< 0* < 1, 



0c + (/if = 0* in Q- 

It is interesting to distinguish in the momentum equation ^ for the cells the contributions due 
to the interaction with the extracellular fluid, represented by the terms (f>c^p and m^c, from the 
remaining ones. In most cases, one can assume that the former are negligible, in terms of order of 
magnitude, with respect to the latter, i.e.. 



(16) 



,\Vp\, |m£c| = o(|V • (0cTc) I, Im^d) in Q, 



so that the main balance is between the intercellular stress and the interaction force with the 
extracellular matrix. Owing to this, Eq. ([6| along with Eq. ( 13 1 provides a particularly simple 
expression for the velocity Vc of the cells: 



(17) 



-V(0cS(0c)). 



By consequence, Eqs. (11), (15 1 give in turn 

1 



(18) 



A. 



V(0cS(0c)) 



K 



-Vp 
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2.2. Mass balance equations. For the system constituted by cells and extracellular fluid, one 
can also write a balance of mass under the assumption that they form a closed mixture, that is 
a mixture which does not exchange matter with the external environment. In more detail, we 
have to take into account that tumor cells, unlike normal cells of the host tissue, are characterized 
by an intense proliferation, which is mainly responsible for the macroscopic growth of the tumor 
cord. In order to focus specifically on this phenomenon, we can reasonably disregard the standard 
reproduction activity of the host cells, so that our mass balance equations read 



(19) 



dt 
dt 



V • ((/.^Vc) = 

V • (0£Vf ) = 



G{^c, c) 
-G(<^c, c) 



in 



(20) 



dt 

d(f>i 
dt 



V-(0,v,) - 
V • = 



In Eqs. (19 1, G{4>c, c) represents a source/sink of cell mass due to reproduction and death of 



tumor cells in connection with the availability of some nutrient c. No contribution of ECM to the 
cell growth is explicitly accounted for, though it is known that the extracellular matrix contains 
growth factors. The reason is that, in the present context, ECM is not supposed to play a primary 
role in the dynamics of the system, hence its inclusion in the function G would result at most in 
constant coefficients related to the volume ratio 0„i (cf. Eq. (14l). 

It is assumed that dead tumor cells are degraded into extracellular fluid, and conversely that 
the latter is consumed whenever a tumor cell duplicates to originate new cells. We incidentally 
pOl) contain the implicit hypothesis that the density of the cells and of the 



notice that Eqs. (19 1 



extracellular fluid is the same, and moreover that it is constant in the whole mixture. 
Denoting by xsi = Xi^(^7 '^)- ^ ^ ^) ^ Q- the indicator function of the set fi, i.e. 



Xn{t, x) 



1 if X G il at time t 
otherwise. 



we can substitute Eq. (17 1 into the first of Eqs. (19 1 to get a single equation for the cell volume 
ratio (t>r'. 



(21) 



dt 



V • 



A. 



!](,/.,)) = G(</)e, c)xn in Q. 



The fiuid volume ratio (f)i is then determined algebraically a posteriori from Eq. (15 1 as 



In addition, summing term by term Eqs. (19 1, (20 1 and recalling again Eq. (15 1 we further 
discover 

V • (0c Vc + (j^i^e) = in Q, 



which, together with Eqs. (17 1, (18), implies 



-V • (KVp) = V 



-V((/.cE(0,)) 



in Q. 



This means that, in view of assumption ( 16 1, we allow for a nonconstant pressure field p that can 
be in turn recovered a posteriori as a by-product of the integration of Eq. (21 1. Finally, once also 
p is known, Eq. ( 18 ) yields the velocity of the extracellular fluid. 



2.3. Cell growth and nutrient diffusion. The function G — G{(l)c, c) appearing in Eqs. (19 1 



(21 1 describes gain (if G > 0) or loss (if G < 0) of cell mass caused by reproduction or death of 
tumor cells, respectively. These are triggered by the presence of some nutrient within the mixture, 
whose concentration per unit volume is denoted by c. As a matter of fact, such a nutrient can be 
mainly identified with the oxygen carried by the blood, which diffuses through the mixture from 
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Figure 3. Logistic cell growth function g{<j)c) (left) and growth regulation func- 
tion r(c) with two thresholds (right). 



the wall of the vessel around which tumor cells grow. It is assumed that oxygen molecules occupy 
no space within the mixture. 

It can be inferred that the basic mechanism by which the available amount of oxygen affects the 
vital functions of the cells consists in the existence of a critical threshold cq > of the concentration 
c. When the latter is above this value, cells are sufficiently fed and can duplicate. Conversely, when 
c falls below cg, cells starve and die. A more sophisticated process may include also a quiescency 
stage, taking place when oxygen concentration is below cq but above a second threshold ci < cq. 
In this state, tumor cells neither reproduce nor die, they simply enter a survival state waiting for 
being reactivated, if c grows above cq, or for definitely dying, in case c falls further below ci. 

/,From the modeling point of view, it is convenient to express the function G in the form 

(22) G(0e, c)=g(0e)r(c), 

so that the new function g carries the information about the specific growth dynamics of the cells, 
while the function F accounts for growth regulation by oxygen concentration. 

We want to point out that a major requirement on the cell volume ratio (j^c is that it should be 
bounded between and 4'.^, , in order to guarantee physical consistency of the solutions returned 



by the model (cf. Eq. (15)). Cell growth mechanisms play undoubtedly a fundamental role in 
this respect, and the choice of g turns out to be crucial in order to generate physically significant 
mathematical models. Therefore, it should rely not only upon physical intuition but also on a 
qualitative knowledge of the mathematical properties of the model itself. Specific features of g 
leading to mathematically well-posed problems will be discussed in detail in Sect. |5] Here we 
simply anticipate that a logistic term (see Fig. |3] left) 

5(0c) = 0c(0* - (f>c) 

is a possible, physically reasonable choice that indeed works. 

Concerning F, it will be clear at the end of Sect. [5] that it mainly influences the distribution of 
tumor cells far from the blood vessel at the steady state. Specifically, it determines the maximum 
size reached by the cord in the transverse x direction. For the moment, we just observe that the 
two scenarios previously discussed may correspond to 

F(c) = 7(c- Co), 7 > 0, 

and to 

r(c) = <( if Cl < c < Co, 70, 71 > 0, 

respectively (see Fig. |3] right). 
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As already mentioned, oxygen reaches cells by diffusion from the vessel wall. Actually, to some 
extent it is also advected through the mixture by the extracellular fluid, but it is known that 
advection is globally quite inefficient with respect to diffusion, due to the relatively slow dynamics 
of the constituents of the mixture compared to the high diffusivity of oxygen molecules (Netti and 
Jain [17'). Therefore we can write for c a standard reaction-diffusion equation of the form 

dc 

(23) — - DAc = -a^cCXii in Q, 

ot 

where Z? > is the oxygen diffusion coefficient within the mixture, and the term at the right-hand 
side models the uptake of nutrient by the sole tumor cells (this is the reason for the indicator 
function xvt) ^ proportional to the product 4>cC via a phenomenological constant a > 0. 

Remark 1. It is worth noting that in the cell growth function G, as well as in the right-hand 



side of the first of Eqs. (20 1, we have not included a term related to the natural death of cells 



(e.g., a term of the form —Scj^c, 5 > 0). Similarly, in Eq. (23 1 neither oxygen absorption by 
normal host cells nor possible chemical degradation of the nutrient (e.g., by a term like t~^c, with 
T > proportional to the characteristic half-life of oxygen molecules) has been taken into account. 
Of course, these effects could be straightforwardly incorporated in the model, but in the present 
context they seem negligible with respect to the main dynamics we are concerned with. 



2.4. Boundary and interface conditions. The system of equations (21 1, (23 1 needs to be 
supplemented by suitable conditions for the two main state variables of the problem, i.e., (j)c and 
c, at the boundary of the domain as well as on the interface S between the regions f2 and f2'^. 

With reference to Fig. [1} we see that the domain Q has four boundaries, whose just one 
is physical, namely the edge Zi representing the blood vessel, while the remaining three serve 
uniquely to confine geometrically the domain of the problem. In particular, I2 can be considered 
as a symmetry boundary, in the sense that one can imagine a symmetric evolution of the system 
for z < 0, while the upper and right boundaries ^3, Z4 should be thought of as 'sufficiently far' 
(ideally, x, z +00) in the host tissue to be unaffected by the dynamics of the growing tumor 
cord. 

In the sequel, the symbol n will denote a generic outward normal unit vector, which has to be 
referred from time to time to the appropriate boundary or interface under consideration. 

At the vessel wall we assume no detachment of both tumor and normal cells, which essentially 
amounts to setting their normal velocity to zero. Moreover, we prescribe the typical oxygen 



concentration found in the blood. Taking Eq. ( 17 1 into account, we have therefore the following 
conditions: 

(24) V((/)cS((/),)) • n = 0, c=Cb on h, 

that in our two-dimensional geometry specialize as dx{4>c'^{<t>c)) = 0, c—cj, for a; = 0. 
On the symmetry boundary we assign canonical symmetry conditions on (/)c and c: 

V(t>c ■ n — Vc • n = on Z2, 

which in our specific case take the form dz4>c = d^c = for z — Q. 

Finally, at the 'far' upper and right boundaries we assume that the host tissue is unperturbed, 
i.e., unstressed, and we set zero flux of oxygen: 

4'c = (f'Qi • n = on Z3, Z4. 

Note in particular that the second of these conditions is interpreted in the present context as 
dxC = on I3, dzC = on Z4. 

According to the classical theory of deformable porous media (see e.g., Preziosi il8j), at the 
interface S between the tumor cord and the host tissue suitable conditions on the continuity of 
the internal cell stress and of the flux of oxygen have to be imposed: 

(25) UcTcuj^O, lVcl-n==0 on S, 

where |-] denotes the jump across the interface S. In particular, given the constitutive relation 
([3]), the condition on the cell stress corresponds to |0cS((^c)l = 0, and, if one assumes that S is 
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a continuous function of 0^, this further reduces to — 0, that is the continuity of the cell 
volume ratio across S. 

In addition, the motion of 5 as a material interface for the cellular matter is regulated by the 
following ordinary differential equation: 

• n = Ve(i, x(t)) • n Vx e 5', 



where Vc is given by Eq. (17), which simply accounts for the fact that, in a Lagrangian framework, 
the points of S move with the normal velocity of the cells themselves. 



3. NONDIMENSIONAL STATEMENT OF THE PROBLEM 

Let i, T, <i>, C denote characteristic values of length, time, cell volume ratio, and nutrient 
concentration respectively, that we choose such that 

L = \/Dt, $ = 0*, C = Cb. 

Let us moreover introduce the non-dimensional variables x, t, <j)c, c defined by the relations 

X = ix, t = Tt, 4>c = (j)*(f)c, C = CfeC, 

the non-dimensional parameters 

^ J, ^0 



and the non-dimensional functions 



1 



9i(t>*4>c), f(c) = rr(cf,c). 



Dropping from now on the tildes from the dimensionless variables and the subscript 'c' from (/),, 
to simplify the notation, the problem originated by Eqs. (21 1, (23 1, along with Eq. (22 1 and the 
boundary and interface conditions (24 1- (25 1, can be formulated in non-dimensional form as 

- V • [0V(0S(0))] - <?(0)r(c)xf2 in Q 



(26) 



dt 



a,(0E(</))) = 0, c=l 
d,(j) = 0, d,c = 
= 00, d^c^O 
(j> = 00, dzc = 
[01 = 0, [Vcl-n = 



in Q 

for a; = 
for z = 
on ^3 
on I4 
on S, 



plus the free boundary condition 
dx(t) 



(27) 



dt 



n= [-V(0I](0))](i, x(t)) -n Vxe5'. 



The choice of 0^ , namely the per cent amount of free space left by the extracellular matrix, as 
reference value for the cell volume ratio implies that at a non-dimensional level we should expect 
< < 1 in Q in view of Eq. (151, as well as < 0o < 1 coherently with the fact that in the 
dimensional formulation of the model 0o is implicitly assumed to be positive and lower than 0* for 
physical consistency reasons. Similarly, since we expect Cb to be the maximum value of nutrient 
concentration found in the system, we consider < c < 1 in Q. 
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t = 325 




t = 650 



t = 900 





0.756 
0.755 
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0.753 
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0.751 
0.75 

0.75S 
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0.754 
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0.75 

0.75S 
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0.75 
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0.751 
0.75 



Figure 4. Contour plot of the cell volume ratio at different times, with values 
increasing from blue to red (range from 0.75 to 0.756). White line indicates the 
position of the interface S. 



4. Numerical results 



Equations (26 1 have been numerically integrated over the reference domain Q = [0, 2.5] x [0, 10] 
via a standard Finite Element discretization, coupled with a level set technique to track the 
evolution of the free boundary S. Specifically, the following functions have been used: 



(28) 

(29) 
(30) 



S(0) 



AP-l 



5(0) = 0(1 - 



r(c) 7(c- Co), 
along with this set of non-dimensional parameters: 

(31) fi^3, 00 = 0.75, 7 = 0.7, cq = 0.8, a = 0.5. 

We refer in particular the reader to Subsect. |5.3|for further motivations on the choice of S. 
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t = 650 0123456789 10 




t = 900 1234567S9 10' 

Figure 5. Contour plot of the nutrient concentration c at the corresponding 
times of the previous Figure |4] with values increasing from blue to red (range 
from 0.7 to 1). White line indicates the position of the interface S. 



Figures |4] [5] show the evolution of the cell volume ratio (j) and the nutrient concentration c, 
respectively, at different non-dimensional times listed on the left. The white line identifies in 
each plot the position of the interface S, thus it visualizes the outer boundary of the tumor cord. 
Initial conditions are set to 0(t = 0) = 0o, c{t = 0) = 1 on the whole domain Q, with ^{t — 0) 
constituted by a quarter of disk of small radius centered at the origin. 

At early times (t = 100) the cord begins to grow keeping a round homogeneous shape, due 
to the high availability of nutrient across the domain. However, the progressive consumption of 
nutrient by the cells makes oxygen concentration fall below the critical threshold cq at the top 
of the cord {t — 325), that is in the farthest zone from the blood vessel. As a consequence, the 
transverse width of the cord starts decreasing {t — 650) and the whole structure takes a more 
elongated shape. At later times {t — 900), a front and a rear regions can be clearly distinguished 
along the cord. The first one is a sort of growing rounded head, where most living tumor cells are 
concentrated because they are globally fed by a sufficient amount of nutrient. The second one is 
a kind of straight tail of nearly constant transverse width, in which the distribution of cells and 
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oxygen depends uniquely on the distance from the blood vessel and has substantially reached a 
steady state. 

The dynamics of the system depicted by the model highlights two main issues to be investigated, 
namely the axial and the transverse growth of the cord. In more detail, the former is concerned 
with the determination of the shape and the growth speed of the head in its motion along the 
blood vessel, possibly in connection with traveling waves describing the evolution of the front part 
of the interface S. The latter is instead aimed at characterizing the typical transverse width of 
the tail, taking advantage of the above consideration that the rear part of the cord is basically 
stationary with an evident special dependence of </) and c on the space variables x and z. In this 
paper, specifically in Sect. [5] below, we address this second problem, leaving the first one for a 
possible forthcoming work. 

We begin from a purely descriptive point of view, observing that in the tail of the cord the cell 
volume ratio attains its maximum value at the vessel wall {x — 0), as it can be expected for in 
that zone cell proliferation is constantly sustained by a direct nutrient supply. Conversely, at the 
interface S and in the corresponding upper region inside the host tissue it decreases toward the 
stress- free value 4>o. Hence the expansion of the tail toward a transverse stationary configuration 
causes tumor and host cells to achieve an equilibrium spatial distribution characterized by a null 
mutual stress as well as a null internal stress of the surrounding tissue. We incidentally notice 
that an analogous unstressed state cannot be reached at the head of the cord, which is never 
stationary for external cells are continuously compressed and pushed away by growing tumor cells 
that penetrate into the host tissue. 

With respect to oxygen distribution, we note that the tail can be ideally divided into two stripes: 
An inner one, extending from the vessel to a certain depth inside the cord, where c > cq, and an 
outer one, reaching the periphery of the cord, in which c < cq. It is straightforward to assimilate 
them to the so-called viable and necrotic regions, respectively, whose existence is experimentally 
confirmed and explicitly described by means of suitable equations in the mathematical model of 
tumor cords by Bertuzzi et al. [5 . However, we want to point out that in the present context the 
formation of these two distinct regions is not postulated a priori as a modeling assumption, but 
is recovered a posteriori as a result of the specific dynamics entailed by the model itself. Oxygen 
concentration, which is maximum at its source, namely the blood vessel, decreases while diffusing 
through the cord due to the absorption by tumor cells. The farther from the vessel cells are the 
lower the quantity of nutrient they are reached by is, so that some of them starve and can no longer 
proliferate. The process continues until in a certain portion of the cord the net proliferation rate is 
zero: Then that portion of cord is in equilibrium, it stops growing transversally and finally reaches 
a steady state. 



5. The stationary problem 



This part of the paper is devoted to the analysis of the mathematical model presented in the 
previous Sects. |2]and[3] In particular, referring to the nondimensional form (26 1 of the equations. 



we consider a steady state in which the cord is supposed to be infinitely long in the longitudinal 
z direction, and to have grown up to a width > in the transverse x direction (see Fig. [6]). In 
such a configuration, the state variables (p and c depend only on x in f2. 

In this geometrical setting, the free boundary S coincides with the line x = w, while the tumor 
cord and the host tissue are represented by the stripe Q < x < w and the half-plane x > w, 



respectively. According to Eqs. (26l, an admissible steady solution (</), c) in Q.'^ satisfies </> — cjjQ, 
Vc — 0, hence we can confine ourselves to the region f2 assuming the continuity conditions cf) ~ (/jq, 
dxC = for X = w. 

The objectives of the present study are twofold: 



(a) As a first step (cf. Subsect. 5.1), to establish existence of physically significant solutions 
with appropriate regularity. In this stage, the width w of the cord is regarded to all 
purposes as a parameter of the model. 
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'far' boundar 



host tissue 




blood vessel 



Figure 6. Geometry used to study the one-dimensional stationary problem. The 
solution (0, c) is assumed to be independent of z, which reduces the domain of 
the equations to the interval [0, w] for the variable x. 



(32) 



(b) As a second step (cf. Subsect. 5.2 1, to solve, via a perturbative approach compatible with 
the results obtained in the previous point, the free boundary problem by identifying a 
steady growth value of the cord width. 

Specifically, we consider the following class of boundary value problem^ 

-dlF{<f>)^w^g{cf>)T{c) ml 

—d'^c + aw'^(j>c — in I 

dx4> = 0, c=l X — Q 

4> = 00, dxC = X = I 

in the unknowns (j), c : I where: 

(i) / = (0, 1) is the rescaled domain after the substitution x — > wx, which reduces the free 
boundary problem to a fixed boundary problem. Notice that w appears now among the 
coefficients of the equations. 

(ii) F ~ F{(l)) is a 'generalized' stress function linked to the cell stress function S = by 
the relation 

(33) F'{<i>) = (/)(0E(0))', 

where the superscript ' stands for derivation with respect to (j). This way it results 

-AF(0) = -V • [</)V(0E(0))] 



in accordance with the first of Eqs. ( |26[ ). 

(iii) g — g{(f>) is the specific cell growth function, which determines how cells proliferate on the 
basis of their current distribution. 

(iv) r = r(c) is the growth regulation function, which expresses promotion or inhibition of 
cellular proliferation as a consequence of the availability of nutrient c. 

(v) a > 0, 00 S (O7 1) are phenomenological parameters related to the consumption rate of 
nutrient by the cells and to the stress-free cell density, respectively. 

In particular, concerning the functions F, g, T we formulate the following assumptions: 
Assumption 1. We assume F £ C^{[0, 1]) with 

F'(0 >0, V^e [0, 1], F'(0 >0, V^e (0, 1], 

and in addition F(0) — 0, F{1) — 1. 



^The condition dx4> = for x = of Problem ( |32[ l is equivalent to the original boundary condition dx{<f>Tj{(f))) = 
of Problem l |26[ l in a sense that will be made precise later in Remark jZj 
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Remark 2. Given the dimensionless cell stress function E, condition F{0) = is readily obtained 



by properly choosing the integration constant in Eq. ( 33 ) . Conversely, we observe that the 



fulfillment of condition F(l) = 1 may simply require a suitable rescaling of E. 



Assumption 2. We assume that g is Lipschitz continuous on [0, 1], with Lipschitz constant 
m2)-9{^i)\<Lg\C2-U VCi, 6 e [0, 1]. 



>0 



Moreover, we require 



giO >Ofor£.e [0, 1], g{0 < otherwise. 



Assumption 3. We assume that T is Lipschitz continuous on [0, 1], with Lipschitz constant 
Lr > 0; 

|ra2)-r(ei)| <irie2-eil, vei,$2e [0, i]. 

From Assumption [l] it follows that F is monotonically increasing and invertible on [0, 1], and 
that its inverse function, henceforth denoted by / for brevity, is in turn increasing and differentiable 
on every compact subset of (0, 1]. We agree to denote by Lf^^ the Lipschitz constant of / on every 
interval of the form [F{e), 1], < e < 1. Notice that the assumptions on F amount mainly to a 
request of ellipticity of the nonlinear equation for (f>, which is however potentially degenerate at 

= 0. 

On the other hand. Assumptions [2] [3] entail the boundedness of g, T on [0, 1], hence there exist 
constants gM, > such that \g{S,)\ < gn, \T{^)\ < Tm for all ^ e [0, 1]. 

For modeling reasons discussed in Sect. |3] we are interested in nonnegative continuous solutions 
(j), c bounded from above by 1. In more detail, in view of the possible degeneracy of the equation 
for (j) in zero, we require that (j) be strictly positive in the domain /, therefore in the sequel it will 
be instrumental to refer to the following sets: 

= {0 e C(/) : £ < (t){x) < 1, Vx e /} , 

U ^ {ce C{i) : < c{x) < 1, Vx e /} , 
where e > is a fixed parameter. Notice that we cannot expect solutions 4> > (pa on /, for they 



would not match the boundary condition of Problem ( 32 1 at x ~ 1 , hence we assume e < 



Remark 3. We anticipate that, as far as the solution to the free boundary problem is concerned, 
the specific value of the parameter e in the range (0, 0o) is irrelevant, because the cell density 
(f> can be shown to satisfy automatically 4> > 4>o ot\ the whole interval / (cf. Theorem [7| . The 
introduction of e is due to the necessity to obtain existence in 14 x J7 of solutions to Problem 



(32 1, among which to look later for possible solutions to the free boundary problem. However, the 
theory we are going to develop will provide as by-product an optimal criterion on whose basis to 
select the value of e. 



Remark 4. The boundary condition for in a; = stated by model ( 26 1 can be rewritten 



according to Eq. ( 33 1 , as 



= a,(0S(0)) = ^^d.cjy. 

Since for </> G 14 one has > e > 0, thus F'{(t>) > in view of Assumption [l] we deduce 
9:,(0S(0)) = O fora; = ^ d^c/) = for x = 0, 



which gives the equivalence of boundary conditions between Problems ( 26 1 and ( 32 ) , at least for 
functions belonging to the class V^. 

In the following, we will consider the main results of our analysis without going into the details 
of their proofs, so as to give an overall picture of the mathematical problem and of its solution. 
The interested reader can find all technical proofs postponed in Sect. |6] We simply point out 
here that the constant Cp, occasionally appearing in several forthcoming formulas, is the Poincare 
constant of the domain /. 
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5.1. Existence and regularity of the solutions. We begin by addressing the question of exis- 



tence and regularity of solutions to Problem (32 1 for fixed positive w. Under suitable assumptions 



on the parameters of the equations, we will identify a family of solutions (0, c) € Vf. x U attached 
to each w ranging in an appropriate set of values. 



Fix then w > 0. The strategy we adopt to solve Problem (32) consists of the following steps: 



(1) We choose any ip £ V^: and study the boundary value problem 



(34) 



- aw'^ipc 



dxc 



in / 



(35) 



establishing existence and uniqueness of a solution c £ U. 

(2) We then choose any a £ U and study the boundary value problem 

-dlF{4)) = w'^g{(t>)V{a) in I 

finding conditions on the parameters that guarantee existence and uniqueness of a solution 

(3) As a consequence of the previous steps, we can define an operator A : ^ Vg that to 
every function g associates the solution cj) G oi Problem ( 35 1 through the solution 



c e C/ of Problem (34). By showing that A satisfies the hypotheses of Schauder Fixed 
Point Theorem, we finally find a function (j) £ V^, and consequently also a function c £ U, 



solving Problem ( 32 1 as desired 



(4) At last, we prove that under the same conditions on the parameters established in the 
previous steps it is possible to control the distance between any solution cj) £ to Problem 
(32) and (j^o, so that a perturbative expansion of (f> about 0o be justifiable in addressing 



next the free boundary problem. 



Steps 1, 2: Well-posedness of Problems (34), (35 1. First, we state the main properties of the 
functions c, (j) s^s they result from any admissible data tp, a on the basis of Problems (34), (35). 



Proposition 1. To every ip £ there exists a unique solution c £ U of Problem (34), which 



is further twice continuously differentiable and monotonically decreasing on I, with a maximum 
value equal to 1 attained for a; = 0. 



Proposition 2. Set 

(36) i3,^p,{s):= 

and define 
(37) 



' Cpgu^M 



P2 — /32(e) CpyJ LgT mLj^^ 



13 := max(/3i(£), /32(e)). 



If Pw < 1, then for every a £ U Problem (35) admits a unique solution (J) £ 



^From the proof of Proposition [2] (cf. Sect. |6|, it turns out that the constants (3i{e), /32(e) 
determine the existence and the uniqueness, respectively, of the solution (p £ to Problem (35). 
Notice that both of them tend to blow when e approaches its limiting values. Indeed, Eqs. ( |36[ ) 
show on the one hand that /3i — > +00 when e — > (/)q", and they may imply on the other hand 
(32 +00 for e ^ 0^ due to the possible degeneracy of Problem (35) for = 0. In fact, observing 



that the Lipschitz constant L/^^ of / on [F{s), 1] can be characterized as 

r 1 1 



(38) 



max f'(£) = max 
«e[F(s),i]' F'(/(e)) 



min F' (rj) 
1] 



we deduce that if F'{0) = then min^gj^ 1] F'{ri) ^ for e ^ 0+, thus L/_e +00. 
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Step 3: Solution to Problem (32 1. We will henceforth always assume I3w < 1, even when not 
explicitly stated. Let us introduce the operator Ai : ^ U that to every ip Q Ve associates the 
unique solution c = Ai {(p) £ U oi Problem ( 34 1 . Analogously, let A2 : U ^ he the operator 
that to such c E U associates the unique solution (fi — A2 (c) e of Problem ( 35 1 with a ~ c. By 
composition, we define the operator A — A2 o Ai : ^ that to every (p £ associates the 
unique solution </> G of the problem 

= «;2g(0)r(c) in/ 

—d'^c + aw'^(pc=0 in I 

dx<j) = 0, c=l X — Q 

(j) = (po, dxC =0 X ^ I. 

Clearly, if A has a fixed point in I4, that is if there exists a function (j) G Vg such that A{(j)) = (j>, 
then the pair (0, c = Ai (</))) e x [/ is a solution to Problem (32 1. Our task is therefore to show 
that A admits such a fixed point. 

For this, it is first useful to know that 

Proposition 3. The operators Ai : U , A2 : U ^ are Lipschitz continuous, and so is 

Moreover, since in applying fixed point techniques a key feature is usually the compactness of 
the involved operator, the following result is particularly welcome. 

Proposition 4. The operator Ai : U is compact, and so is A : Vg V^. 



Thanks to Propositions |3j |4| we can apply Schauder Fixed Point Theorem to solve Problem ( 32 1 . 
For the sake of completeness, we explicitly recall in the statement of the theorem all hypotheses 
that bring to the result. 

Theorem 5. Let Assumptions Q] [?| hold with Pw < 1, where (3 — (3{e) is the constant defined 
by Eq. (37 1. Then there exists a solution (</>, c) £ Vi; x U to Problem ( |32[ ) satisfying the a priori 
estimate 

||1 ~ c||oo < aw'^\\(j)\\L^i)- 

Theorem [5] gives existence but not uniqueness of solutions {(f>, c) € x U to Problem (32 1. 
Observe however that, owing to Proposition [l] to any fixed point (j) E Vg of the operator A there 
corresponds a unique c = Ai{4>) G U such that the pair (0, c) G Ve x U solves Problem (32). 

Step 4- Controlling the distance between (j) and 0o. For any solution (cj), c) G x U to Problem 
(32 1, we are interested now in controlling the distance between </> and the constant (j)Q, in view of 
the application of a perturbative technique to solve the free boundary problem. We find out that 
the quantity P2W plays a fundamental role, in fact we have: 

Theorem 6. Let (p G V^: be any solution to Problem ( p2| ). Then 

u~Mo.<^{f32wr. 

^From this result we infer that, since [32W < j3w, for even moderately small values of (3w (with 
respect to 1) the solution </> becomes close to (jj^ provided qm (maximum cell proliferation) and 
Lg (maximum proliferation rate) satisfy gu/Lg = 0{Cp), Cp ~ 0.64 (cf. Sect. [6]). 

5.2. Solution to the free boundary problem. So far we have obtained existence and regularity 
of solutions to Problem ( 32 1 under a suitable condition on the admissible steady values of the cord 
width w. However, we are mainly interested in finding stationary solutions satisfying in addition 
the free boundary condition (27 1, that in the present setting rewrites as 

d^{(j)Y.{<j))) ^ forx^l 

or, on the basis of the same considerations proposed in Remark [4] more explicitly as 

(39) d,j:(j) = for a; = 1. 

To this end, we formulate additional assumptions on F: 
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Assumption 4. In addition to Assumption^ we assume that T is nondecreasing on [0, 1]; 

r(6)<r(6), va,6e[0, 1], a<6 

with r(o) < < r(i). 

Notice that Assumptions [s] and |4] imply the existence of a point cq G (0, 1) such that r(co) = 0, 
with moreover r(^) < for ^ € [0, co), r(^) > for ^ e (cq, 1]. 

We state first of all a crucial result that clarifies the role of the parameter e. 



Theorem 7. Let Assumption^hold and assume that a solution 
(32 I, such that (f> fulfills the free boundary condition (39 I. Then i 



t>, c) Ve y< U exists to Problem 
(x) > (f>Q for all X £ I . 

In view of Theorem]?] the specific value of e S (0, (f>o) is not relevant to solve the free boundary 
problem, indeed any possible solution, provided it exists, is 'naturally' bounded from below by (f>o. 
Considering that Theorem |5] guarantees existence of solutions if (3w < 1 and that, on the other 
hand, it results /3 +oo when both e ^ 0+ and e ^ (p^ (cf. Proposition |2]) , we deduce that 
the optimal criterion for the choice of s aims at minimizing /3(e), so as to get the widest possible 
range of admissible values of w. 

Let us look for a first order expansion of (/) in a neighborhood of the constant stress-free value 



(40) 



where > is a 'small' parameter to be identified from Problem (32 1. Correspondingly, also c is 
expanded with respect to in a similar way: 



(41) 

Inserting Eqs. 



c{x) 



(x) 



vc 



(1), 



x). 



(40), (411 into the second of Eqs. (32 1, and retaining only the terms up to the 



= 
= 1 
= 



in / 

X = 
X — 1, 



zeroth order in yields the following problem for c' 
(42) \ c(0) 

whose solution reads explicitly 

cl°H^) = 



(43) v-^; - W / ^ ^ 

cosh [w^JacpQj 

where we have emphasized the dependence of c'"' on w as a parameter. It is an easy task to check 
that Cw'^ G U, indeed Problem (42 1 falls into Proposition [l] with the special choice ip = 4>o & Ve- 
in order to get an equation satisfied by (t>^'^\ and to detect at the same time the parameter 
v, we proceed in a similar fashion, introducing expressions (|40|), (41 1 into the first of Eqs. (32 1. 



x)) 



Assuming formal differentiability of g, T we find, after some standard algebra. 

As r is not identically zero, another term of zeroth order in should be found in this equation, 
besides the one at the right-hand side. Estimating the order of magnitude of T{cw'') via Tm and 
looking then at the left-hand side we discover v = 0{w'^T m / F' {(f'o))- Moreover, recalling the 



expression (38 1 for the Lipschitz constant L/.^ of / on [F{e), 1], we see that, at least for suitable 
functions F, we may have F'{(j)o) = 0{LJ^^), hence 

(44) ly = Oiw^TMLf,,) = Oi{p2wf). 

As anticipated above, the parameter v must be 'small' in order for the perturbative expansion 
to be meaningful. Equation ( |44| shows that this requirement corresponds in essence to the same 
condition prescribed by Theorem [s] for the existence of solutions to Problem (32 1. Therefore, the 
existence theory and the perturbative approach are consistent with one another, as they impose 
the same conditions on the parameters w of the model. In other words, for a certain value of 
the quantity l3w, they either apply or fail both. 
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Retaining only the terms up to the zeroth order in v, the problem for the perturbation 0^^^ 
reads: 

= C5(0o)r(cl°)) ml 
= x^l. 



(45) 



where C — 0{Tj^j) is an appropriate positive constant. In particular, boundary conditions for 
0^^) are deduced from those prescribed on (j) taking Eq. (40 1 into account. 



Unlike the previous case for cff ^ , Problem ( 45 ) cannot be solved in a closed form for a generic 



function T. However, classical theory for linear elliptic equations guarantees existence and unique- 
ness of a solution cf)^^'' G Hq which is actually continuous due to Sobolev embedding theorems 

for / CM. Furthermore, since Cw'^ G U and T is continuous on [0, 1], we even get d'^4>^^^ £ 



hence (j)^^^ is a classical solution to Problem (45 1. This allows us to evaluate its derivative for 



X & I hy integrating once the differential equation in ( 45 1 



-dA^'\x) + dA^^'\Q) = Cg(0o) j r(4°)) d^, 



whence, owing to the boundary conditions for x — Q, x — 1, 



(46) 



T{c^^'^)dx^Q. 



Solving the (approximate) free boundary problem corresponds therefore to finding the roots w > 
ofEq. (liel). 



After some preliminary considerations about the dependence of Cw^ on w 



Proposition 8. Let c 



(0) 



(zU he the solution to Problem (42 I. Then: 

.(O)t 



JO) 



(x) 



(i) The mapping w i— > Cw (x) for w G M+, x £ I is continuous and nonincreasing. 

(ii) For w — > +oo the functions cL?' tend pointwise to the function 

' 1 ifx = 
ifx^O. 

(iii) For all cq € (0, 1) there exists > such that if w > w^, the equation 

in the unknown x has exactly one solution x^ G I- 

(iv) For any fixed cq G (0, 1) one has — *■ when w +oo. 

we have all we need to solve our free boundary problem. The following result should be compared 
to the analogous one proved by Bueno et al. |10j . 

Theorem 9. Let Assumptions\3\\4\hold. Then there exists a positive solution Wq to Eq. (46). 



Notice that Theorem [9] does not guarantee that the solution wq actually satisfies Pwq < 1. In 
practice, this condition has to be checked a posteriori, after solving Eq. (46 1 explicitly. Concerning 



this, we observe that, once the function F has been specified, Eq. (46 1 often specializes in an 



algebraic equation in the unknown w, that can be solved with arbitrary precision, possibly with 
the aid of suitable numerical techniques. After getting the solution wq claimed by Theorem [9| two 
situations may arise: Condition f3wQ < 1 is either satisfied, hence, in view of the reasonings above, 
one can take wq as a good approximation of the solution to the free boundary problem, or it is 
violated, in which case one should reject wq and conclude that, for the specific model parameters 
at hand, no information on the existence of a solution satisfying the free boundary condition can 
be obtained from the present theory. 

We show an example of application of this procedure in the next Sect. |5.3| 
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Figure 7. Cell stress function of Eq. (47 1 for ^ = 1, 2, 3. 



5.3. A noticeable example. In order to illustrate how the previous theory applies to a specific 



model, we consider for /i > 1 the following class of intercellular stress functions E (cf. Eq. (28l) 



(47) 



■log- 



if ^ = 1 



if ^ > 1, 



whose behavior for several /j, is illustrated in Fig. [7] Notice that for /i > 2 the function E 
grows unboundedly with </>, and is thus compatible with the expected behavior of an intercellular 
pressure-like stress, while for 1 < /i < 2 it is bounded from above, and if ^ 2 it even tends to 
zero for (j) — + +oo. From the modeling point of view, the meaningful cases correspond to > 2. 
However, the theory developed in the previous sections also covers the range of values 1 < /i < 2, 
which may be of some theoretical interest due to the particular form of the equations they originate 
(for instance, for fJ. = I the differential part of the equation for (j) is linear). 



According to Eq. ( 33 1 , the corresponding generalized stress function F is 



F(0) = r, 

which gives rise to the stationary porous medium equation for the cell volume ratio </>. It is 
straightforward to check that such an F complies with Assumption [T] for all /x > 1. Moreover, we 
choose for g and T the expressions given by Eqs. (29 1, (pO]): 



5(0) = 0(1-0), r(c) = 7(c-co), 

where < cq < 1 and 7 > are constant, whence also Assumptions [2] [3j [4] are satisfied. 
The inverse function / of on [0, 1] is 

/(0) = </0. 

Notice that / is actually Lipschitz on every compact subset of (0, 1] for all /i > 1, but for /i > 1 
it is not Lipschitz on the whole interval [0, 1] due to the vertical tangent at = 0. However, we 
know that the interesting parameter is its Lipschitz constant Lf^^ on [F(s), 1], which exists for 
any M > 1 and equals 

1 1 



max f'(f) 



max 
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Figure 8. Numerical solution (FEM) to Problem (49 1 fulfilling the free boundary 



condition (39 1. Cell density (f) on the left, nutrient concentration c on the right. 



Moreover, 



5M= maKjff(ai = 4, 



L 



9 - ™ax^ \g\C}\ 



Ce[o, 1] 



Tm = max \T{Ci\ = 7max (cq, 1 - cq), 
«e[o,il 



whence we compute 
/3i(£) 



/7max(co, 1 - cq) 



277 1 



en 



2 jma,x{cQ, 1 — cq) 
77 y /xe^'-i 



For all ui > such that Pw < 1, that is 
(48) w < min 



1 1 

Theorem [5] guarantees the existence of a solution (</>, c) g x U" to the problem 

—dl4>'^ — 7W^0(1 — (/))(c — Co) in / 
—d^c+aw'^(j)c—0 in I 

dx4i — Q,c—l X — Q 

4) = 00, =0 x^l. 



(49) 



Let us now look for a solution to the free boundary problem ( 39 1 via the small parameter 



approximation discussed in Subsect. |5.2[ The specific form of the function F we are using allows 



( 43 1 into the latter gives 



for an easy estimate of the solutions to Eq. (|46|. Inserting the expression of c^u^ given by Eq, 



(50) 



tanh 



{wyja4)i^ = w\/a0oco, 



thus, invoking McLaurin expansion of the hyperbolic tangent and considering that tanh(a;) < 1 
for all X e M, 



(51) 



/ 3(l-co) 1 

T < Wo < 



For the sake of definiteness, let us fix the same parameters already used in Eq. (31 1, then let 
us determine the optimal e by solving the min-max problem 

min max Bi (e) . 

e6(O,0o)*=1.2 
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Figure 9. Relative errors on (left) and c (right) as they result from comparison 



between the approximate solution, analytically computed via Eqs. (43 1, (45 1, and 



the numerical solution (FEM) to the exact equations (Problem (|49|)). Notice that 
in both cases one has E = 0{1Q~^). 



This gives e ~ 0.5 and correspondingly /3i ~ [32 



0.55, whence it can be easily computed that the 
Solving Eq. ( 50 1 via an iterative numerical 



estimate (51) agrees with the previous bound (48 1 
procedure yields wq ~ 1.45. Since j3wQ ~ 0.80, the small parameter assumption and the related 
perturbative approach are consistent with the problem at hand, as it is definitely confirmed by 
the value of w returned by the Finite Element solution (see Fig. |4] [8 
w ~ 1.44. The behaviors of the relative errors on (f) and c (see Figure 



of the exact problem, i.e.. 



E[4>]{x) 



1 



0o + t.0(i)(a;) 



E[c]{x) 



JO) 



1 



c{x) 

where (j)^^"^ is analytically computed for this particular case from Eq. ( [45| using C — TJj, and 
where we have set v — {/32w)'^ , further prove the good quality of the small parameter approximation 
with respect to the exact solution. 



6. Technical proofs 

Here we collect, without specific comments, the proofs of all results stated in the previous 
Sect. [5j Before going into technical details, however, it is convenient to fix some basic facts and 
notations that we will extensively use in the sequel. 

We denote by Hq q{I) {Hq i{I), respectively) the closed Hnear subspace of the Sobolev space 
H^{I) = W^'^{I) consisting of all functions u G H^{I) whose trace vanishes at a; = (x = 1, 
respectively). We recall that Poincarc inequality holds true for functions u E Hqq{I), or u G 

where the Poincare constant Cp equals ^. As a consequence, the quantity 

— \\dxU\\L^I) 

defines a norm in both Hqq{I) and Hq ^{I) equivalent to the usual i/^-norm. Notice that in 
referring to such a norm we use generically the short subscript H(I) for either of the spaces, the 
meaning being recoverable each time from the context. 

When dealing with continuous functions u £ C{I), we indicate by ||w||oo their oo-norm over /: 

||w||oo = max |u(a;)|. 

Finally, we use the symbols for the positive and negative part of a function u: 

=max(u, 0), =max(— u, 0). 
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In view a theorem due to G. Stampacchia, it is known that u e H^{I) imphes w^, u G H^(I) 
with moreover 

I \da;U ifu>0 _ \-dj;U ifu<0 

I otherwise, I otherwise. 

Let us start by a mainly technical Lemma, which will be sometimes referenced in the forthcoming 
proofs. 

Lemma 0. Let h, k £ C{I) with moreover h{x) > for all x £ /. If u £ C{I) solves the equation 

—d^u + h{x)u — k{x) in I 
with either of the following sets of boundary conditions: 

(i) u = for X — 0; dxU = for x — I, or; 

(ii) dxU — for x = 0; u = for x = 1, 
then 

Hx)\<\\k\\^, Vxe/. 

Proof. Using standard techniques for linear elliptic problems, it can be shown that the proposed 
equation admits a unique solution u G q (/) in case (i) , and a unique solution u € Hq ^ (!) in 
case (ii). In both cases, u G C{I) due to Sobolev embedding theorems for / CM. Moreover, 
Lax-Milgram Theorem entails the following a priori estimate: 

\\u\\h{i) < \\k\\L--{i), 

which, considering that 

1 1 

2 



= / \kix)\'^dx< I |^max|fc(a;)| ) dx ^ \\k\ 



oo 



and, from Morrey inequality, that ||it||oo < yields 

||"||oo < ll^lloo 

and thus the thesis. □ 



Proposition 1. To every tp G there exists a unique solution c G U of Problem (34 1, which 
is further twice continuously differentiate and monotonically decreasing on I, with a maximum 
value equal to 1 attained for x = 0. 



Proof. (1) By means of the substitution c = c — 1, Problem (34 1 is converted into the following 
linear elliptic problem with homogeneous boundary conditions: 

-dxC + aw^ipc — —aw'^ip in I 

c = a; 

dxC — X = 1. 

By putting it in weak form as 

find c G H^q{I) such that 

a{c,v)^b{v), yvGHloil), 

where the bilinear form a : Hq g(/) x Hq M. and the linear form b : Hq q(/) M. are 

defined by 

11 1 
a{u, v) — / dxudxvdx + aw^ / (puvdx, b{v) = —avo^ f ipv dx 



(52) 



respectively, classical theory for linear elliptic PDEs can be applied up to checking that a 
is continuous and coercive on Hp g(/) x //q g(/) and that b is continuous on Hq q{I). Owing 
to Lax-Milgram Theorem, we get existence and uniqueness of a solution c G Hqq{I) to 
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Problem (52 1, and by consequence of a solution c e H^{I) to Problem (34 1. Moreover, for 
/CM Sobolev embedding theorems imply H^{I) C C(J), thus c is in fact continuous on 
/. 

(2) We check now that c > on /. For this, observe first that £ Hqq{I), then multiply 
Eq. (34) by and integrate by parts over /: 



I if (7) 



aw"^ / (p{c ) dx = 



whence 



||c \\h{i) — —aw / ip{c ) dx < 



We deduce |lc^||/f(/) = 0, which yields = almost everywhere in / and, by continuity, 
c{x) > for all x £ I. 

(3) Finally, we show that c < 1 on /. Since c e C(/), the second order derivative d^c exists as 
a distribution. But from Eq. (34 1 we discover that it actually coincides with a continuous 
function: 

d^c = aw'^tpc e C(/), 



so that c turns out to be a classical solution of Problem (34 1. This allows us to compute, 
for every x € I, 

1 



id,j:c){l) - {dxc){x) = aw^ J ipcd^, 

X 

whence, using the boundary condition at x = 1, 

1 

id^c){x) = -aw^ (pcd^<0, Va; G /. 



Therefore we have that c is nonincreasing on /, and consequently we deduce 



maxc(a;) = c(0) = 1, 



which completes the proof. 
Proposition 2. Set 



□ 



/3i=/3i(e) 



CpgM^M 



and define 



F{M - F{e) 
(3 = l3is) max(/3i(e), /32(e)) 



/32 = /32(e) -.^Cp^LgTMLf,, 



If Pw < 1, then for every a £ U Problem (35) admits a unique solution (j) £V^. 



Proof. (1) First we rewrite Problem ( |35[ ) making the substitution 
(53) u^F{<P)-F{c^^), 

which in turn implies (j) = f{u + uq) for uq :— F{<j)^) e (0, 1): 
—d^u ~ w'^g{f{u + u^))T{a) in/ 
dxU = a: = 

u = Q X ~ \. 



(54) 



Notice that for </> S the boundary conditions at a; = of Problems (35 1 and (54) agree, 
since d^u = F'{4>)dx(j) with F'^cj)) > 0. 
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(2) Now we prove that for all it; > Problem ( 35 1 admits continuous solutions ranging in 
[0, 1]. In view of Eq. (53 1 and of the continuity and monotonicity properties of F, this 



(55) 



amounts to showing that Problem (54 1 possesses solutions u e C(/) bounded between 
—Uq and 1 — itQ. 

(a) Let us introduce the function 

~g{u)^[g{f{u + uo))] + 
and consider then the auxiliary problem 

-d^u = 'w^g{u)T{(7) in / 
dxU = X = 

u ~ X — I. 

The function g{u) being continuous, if P{u) denotes any of its antiderivatives on M 
then Problem (551 can be viewed as the Euler-Lagrange equation associated to the 
energy functional E : Hq ]^(/) R, 



E{u) 



1 1 
{d^ufdx-w^ J P{u)T{a) dx. 





Due to the assumptions on the sign of g and the monotonicity of /, it results g{u) = 
for u < —uq or u > 1 — uq, hence P{u) is bounded: There exists a constant C > 
such that |P(u)| < C for all u e M. This entails first that E is finite on H^ i{I): 



E{u) < 



1. 



and moreover that it is coercive: 



+ w^CTm < +00, 



E{u) > 



> ^\\u\ 



\H(I) 



2 

\h{I) 



\P{u)\-\T{a)\dx 



w^CT 



M 



> 



where C" > is a new suitable constant. Thus any minimizing sequence is bounded 
in Hq and has therefore weakly convergent subsequences. If we can show that E 
is sequentially weakly lower semicontinuous on Hq we will deduce the existence 
of a minimizer u € Hq that is a solution to the auxiliary problem (55 1. For this, 
notice that E belongs to the class of functionals of the form 



E{u) 



L{dxU, u, x) dx, 



where L : x / M is the function given by 

w^P{z)T{a{x)). 



Hp, X) = ^ 



Clearly, L is convex in p for each z e M, a; G /, and moreover it is bounded from below 
as it results L{p, z, x) > —liP'CV m for all (p, z, x) £ x /. These two features are 
sufficient in order for E to be sequentially weakly lower semicontinuous on H^{I) (see 
Evans [15 , Chapter 8, p. 446 for further details), hence also on H^ ^{I) as desired. 



We conclude that the auxiliary problem (551 has solutions u E Hq which, as 
a consequence of Sobolev embedding theorems for / CM, are actually continuous 
functions u G C{I). 
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(b) We claim that any solution u e Hq ^{I) to Problem (55 1 satisfies the inequalities 
—uq < u{x) < 1 — uo for all x £ I. To see this, let us multiply Eq. (55) by 
(u + uq) G Hq and integrate by parts over /: 

1 

dxudx{u + uo) dx = J g{u){u + uo) T{a)dx. 



Since g{u) = for u < —uq while (u + Uq) = for u > — ug, this relation implies 

||(w+Mo)"llff(/) = 0, 

hence, by continuity, u(x) > — uq for all x Cz I. 

Analogously, multiplying Eq. (55) by (1 — wq — w) G Hq and integrating by 
parts over / gives 

1 

dj:U dx{l ~ Uq — u) dx = J g{u){l — UQ~u) T{a)dx. 



But g{u) — for u > 1 — Uq while (1 — — m) =0 for u < 1 — uq, which entails 

11(1 - Uq - uy\\H(I) = 

and finally u(x) < 1 — uq for all x Cz I. 

(c) We observe now that for u e [~uq, 1 — uq] one has g{u) — g{f{u + uo))i thus any 
solution to the auxiliary problem is also a solution to Problem (54). Therefore, the 
latter admits solutions u e C(/) bounded between — uq and 1 — uq, which, via Eq. 
(53), originate solutions ((> G C{I) to Problem (35 1 such that < < 1 on /. 

(3) Finally, we show that it is possible to get solutions (/) S to Problem ( 35 1 provided w is 
sufficiently small. For this, let g C(/), < < 1, be a solution to Problem (35). Since 
F{(j)Q) is constant, we can rewrite the equation for cj) in the equivalent form 

-dl{F{4^)~F{M)^y^^gm{<^)- 

Observe now that F{(l)) — F{<j)o) £ Hq with moreover 
dx {F{^) - F{^q)) = F'{cf>)dxcf> = 

at X — 0. Thus, multiplying both sides of the above equation by — F{(l)o) and 

integrating by parts over /, we discover 

\\F{(b) - F(0o)||h(7) < w^CpgMTM, 

where we have used the fact that, for (j> G [0, 1], the function g is bounded by gM- In view 
of Morrey inequality we further deduce 

\\F{q}) - FiMWoo < w^CpgAiTM, 

whence 

F(0o) - w^CpgMTM < F{4>) < F{cf,o) + w^CpguTM 

on /. Condition (p > e is equivalent to > F{e) due to the monotonicity of F. 

Imposing 

F{e) < F(0o) - w^CpguTM, 
we obtain G Vg provided Piw < 1, where 



/3i=/3i(e) 



CpgM^M 



^ Fi^o)^F{s)- 

Notice that Pi is well defined, since e < 0o implies F{e) < F{(j)Q). 
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(4) Regarding uniqueness, suppose, under the hypothesis (3iw < 1, that i/ii, 02 G two 
solutions to Problem ( 35 1 . Subtracting the respective equations we find 

-9^(F(02) - = i«2[g(02) - 5(0i)]r(a), 

where ^"(02) - F{(l)i) e i?o,i(/) is such that 

for a: = 0. Multiplying both sides by F{(j)2) — F{<j)i) and integrating by parts over I we 
get 

1 

2 _ „,,2 



||i^(02) - F{<j,,)\\ij^r) j [9{<i>2) - 5(0i)](J^(02) - F{cj,{))T{a) dx 



1 

< w^LgVM j \4>2-4>i\- \F{,4>2) - F(0i)| dx 





1 

,,2 





< w^LgVMLf.,,Cl\\F{(p2) - F{^,)\\],^,), 
where we have used the Lipschitz continuity of g on [0, 1]. Setting 

h = /32(e) = Cp^LgVMLf.,, 

we deduce therefore 

{\-pW)m<ly2)-F{(ly^)\\l^,)<Q, 

thus if w satisfies the further constraint /32U' < 1 we conclude 

||F(02)-F(0i)||h(/) -0, 

that is F{(j)2) = F{(j)i) for all x ^ I, and the uniqueness follows from the bijectivity of F. 
(5) In view of the previous results and given the definition of /3, we conclude that if f3w < 1 
then there exists a unique solution € 14 to Problem (35 1 as desired. □ 

Proposition 3. The operators Ai : U , A2 : U ^ are Lipschitz continuous, and so is 

Proof. (1) Let us begin by considering Ai. Given ipi, ip2 & let ci — Ai{ipi), C2 = Ai{ip2) 
be the corresponding solutions in U to Problem (34 1 and set u = C2 — ci £ C{I). Then u 
solves the problem 

—d'^u + aw'^ip2U = —aw'^ci{ip2 ~ ^pi) in/ 
u = x^O 
dxU = a; = 1. 

Setting 

h{x) = aw'^(p2{x), g{x) = -aw'^ci{x){ip2{x) - >fi{x)), 
we have h, g d C{I), h> on I, and moreover 

II.9II00 < aw^\\'P2 - <y5l||oo 

because Ci G U satisfies |ci| < 1 on J. ^^From Lemma [o] we deduce then 

\\Ai{ip2) - Ai{ipi)\\ao = \\C2 - Cilloo = ||w||oo < llglloo < aw'^\\ip2 ~ (^lHoo, 

whence the thesis follows. 
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(2) Concerning A2, let CTi, (T2 G ?7 and let (f>i = A2{ai), 02 = ^2 ('''2) be the solutions in to 



Problem ( 35 1 , respectively. Subtracting the corresponding equations we find 

-dl{F{cf>2) - F{cl,,)) = [5(<^2)r(a2) - 5(0i)r(ai)] 

= [5(</)2) - g{^i)] r(a2) + w^g{(f>i)[T{a2) ~ T{ai)]. 

After observing that F{(j)2) - F(0i) e i?o,i(^), d^{F{4>2) - F{(t)i)) = for x = 0, we 
multiply both sides by F{(j)2) — and integrate by parts over /, like in the proof of 

Proposition |2] to obtain 

1 

\\F{cI,2)-F{cI>,)\\1^j^<w^Tm J \g{<l>2)-9{4>i)\-\F{cl>2)-F{cl,i)\dx 



1 

+ w' [ |5(0i)| • |r(a2) - r(ai)| • \F{(b2) - F{q^,)\dx 





<Plw^F{^2)^F{^,)\\j,^j^ 

+ w'^gMLrCp\\(J2 - o-i||i2(/)|li^((/)2) - F{(j)i)\\H{i)- 

Notice that in handling the second integral at the right-hand side we have used Cauchy- 
Schwarz inequality in L^{I). Thus 

(l-/32'u;2) 11^^(02) < 

w'^gMLrCp\\(T2 - cri|lL2(/)||F(</i2) - F{(j>i)\\H{i), 
whence, recalling that P2W < Pw < 1 and dividing by \\F((j)2) ~ -^(01)11^(7) if 0i 4>2, 

||F(02) - F{(t)i)\\H(I) < '^,^'"m'^2^ \\^^ - Cri||i2(j). 

Furthermore, by Morrey inequality and using ui, a2 ^ U (1 C{I) we get 
Hf(<^2)-^(0i)||oo< "'%t? lk2-a,|U. 

i P2^ 

But 

\(f>2 im02) - m0l)l < Lf,e\Fi(t>2) - F{^i)l 

so that we finally have 

\\A2{<J2) - A2{a,)\\^ = U2 - Moo < ^'g^^-^rCp-L^.e n^^ _ ^^11^ 

1 — P2?i' 

as desired. 

(3) Lipschitz continuity of A follows immediately by composition. □ 

Proposition 4. The operator Ai : U is compact, and so is A : V^. 

Proof. (1) Continuity of Ai is implied by Proposition [3] hence in order to have compactness 
we need to prove that Ai{Ve) is relatively compact in U. We do it by showing that for 
any {ipn} C the sequence {c„} — {Ai(ipn)} C U contains a convergent subsequence. 



Due to ||(^„||oo, lic„||oo < 1, Eq- ( |34| yields ||(3^c„||oo < cxu)-^ for all n, hence, using the 
boundary condition at x = 1, 



1 

|2, 



9.c„(x)| < / \d'^cnmd^<aw' 



MULTIPHASE MODELING OF TUMOR CORDS 



29 



In view of these estimates, we conclude that the sequence {c„} is uniformly bounded and 
equicontinuous. Owing to Ascoli-Arzela compactness criterion, we can therefore extract a 
subsequence {c„^} converging uniformly on / to some c € AiiVf,): 

lim ||c„^ - c||oo = 0, 

k — *oo 

whence the compactness of Ai follows. 
(2) To show the compactness of A we rely on that A = A2 o Ai. Continuity of A is implied 
by Proposition [3] Moreover, if {ipn} is a sequence in Vg and we let c„ — Ai{tpn), 4>n — 
A2{cn) = A{(pn), the compactness of Ai previously proved allows us to assume, passing 
to a subsequence if necessary, that {c„} converges in U. But then the continuity of A2 
implies that also {4>n} converges in V^, and we have the thesis. □ 

Theorem 5. Let Assumptions^ [5| /lo/d with j3w < 1, where (3 = (3{e) is the constant defined 
by Eq. (37 1. Then there exists a solution ((/>, c) £ x U to Problem ( |32[ ) satisfying the a priori 
estimate 

||1 - c||oo < aw'^\\4>\\L2(jy 

Proof. We know from Proposition [4] that A : is compact. We claim now that is convex 

and closed in C{I). 

(1) Convexity: Let u, v £ V^, X G [0, 1], and define 

sx{x) — Xu{x) + (1 — X)v{x), X G I. 

Since u, v are continuous so is s\, and moreover: 

(i) sx{x) > Ae + (1 - A)£ = £ 

(ii) sxlx) ^ A + 1 - A = 1 

for all x G I, hence s\ G for all < A < 1. 

(2) Closure: Let {u„} C K be a sequence converging to some u g C(/), that is — u|loc — > 
when n — !■ 00. Then converges pointwise to u on /, hence we must have 

u{x) — e — lim {un{x) — e) > 0, 

n — >oo 

and analogously 

1 — u{x) = lim (1 — Un{x)) > 0. 

n — >QO 

We conclude therefore e < u{x) < 1 for all x G I, that is u G V^. 

According to Schauder Fixed Point Theorem, A has a fixed point <j) G Ve- Setting c = Ai{(j)), 
we deduce that the pair {(f>, c) G x U solves Problem (32 1. 

The a priori estimate on 0, c readily follows from Lax-Milgram Theorem applied to Problem [T| 
with If ^ (p. □ 

Theorem 6. Let (j) G be any solution to Problem ([32|). Then 



Proof. ^From Proposition [2] we know 

\\F{(j,) - F(0o)lk(7) < w'^CpguTu. 
Rearranging the coefficients according to the definition of ^2 given by Eq. ( 36 1 we get 

LgLf,,\\F{4>)-F{M\\H(,i) < ^(/32w;)^ 

Op 

whence the thesis follows using Morrey inequality at the left-hand side and then considering that 

\^ix) ~ < LfjFi^ix)) - F{^o)\ 
for all a; G / because (/)> e on I. □ 

Theorem 7. Let Assumption^hold and assume that a solution ((/>, c) G x U exists to Problem 
(32 1, such that (p fulfills the free boundary condition (39 1. Then 4>{x) > (j)o for all x G I. 
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Proof. (1) Let us integrate once the differential equation for (f> on /. Owing to the boundary 
condition dx(t){0) — and to the free boundary condition dx(f>{l) = we find 

1 

g((/))r(c) dx = 0. 

Since g{(t>) > is continuous and cannot identically vanish on / (basically because the 
constant ip = 1 does not solve Problem ( 32 ) , as it does not match the boundary condition 
0(f) — (po < 1) , we deduce that either T{c) = on / or T{c) changes sign on /. Observing 
that, due to the properties of c (cf. Proposition [IJ and of F (cf. Assumption |4]) , the 
function x '—^ r(c(x)) is continuous and nonincreasing on / with r(c(0)) = r(l) > 0, we 
conclude that there must exist x E I such that r(c) > in [0, x) and T{c) < in [x, 1]. 

(2) From the equation 

-dlF{<j>) = w^g{<P)T{c) in / 

we discover that 

' < for x e [0, x) 
> for X e [x, 1], 

hence that F((f)) is concave in the interval [0, x] and convex in the interval {x, 1]. Since 
dxF{(j)) = F'{(j))dx(i), we further obtain (5^F(<?i))(0) = 0, whence dxF{(j)) < in [0, x), 
and {dxF{(t))){\) — 0, thus dxF{(j)) < also in [x, 1]. In conclusion, we have dxF{(f) < 
on the whole /. 

(3) Assume now by contradiction that (F((/)))(a;) < F((/)o) for some x E [0, 1). Then 

1 

(F(0))(1) = (F(0))(x) + / dxF{c^) < FiM, 



dim 



which however is incompatible with the boundary condition 4>{1) = 0o of Problem (32 1. 
Hence F(0) > F(0o) on /. 
(4) Finally, we see that it is impossible to have (l){x) < (po at any x € I, for this would imply 
(F(0))(a;) < F{(f>o) which contradicts the previous result. Therefore > 0o on / and we 
have the thesis. □ 



Proposition 8. Let c^^ £ U be the solution to Problem ( |42| . Then: 

(i) The mapping w Cw\x) for w G M+, x G I is continuous and nonincreasing. 

(ii) For w +oo the functions Cw'^ tend pointwise to the function 



1 ifx = 
ifx^O. 

(iii) For all cq G (0, 1) there exists > such that if w > w^, the equation 

in the unknown x has exactly one solution x^ G /. 

(iv) For any fixed cq G (0, 1) one has x^j — > when w +oo. 



Proof. Some of these properties are more easily proved by referring to Problem ( 42 ) rather than 
to its explicit solution ( 43 1 . 

(1) To show continuity of the mapping w Cw\x), let us fix wi, W2 > and set u = 
Cw} — Cw} G C(/). Then u solves the problem 

—dxU + aw24>QU = —a{w2 — wI)(I)qCwI in / 
(56) < u = x^O 

dxU = X = 1 
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SO that, owing to Lemma [Oj we have 



(0) - cl°)|U = hiloo < aM^j - u;?||!<)|U < a^wi - wl\. 



w-i I loo 



Consequently Cwl{x) —> Cw}{x) for all a; G / as W2 — > wi, and continuity follows. 

(2) We address now monotonicity ofw Cw\x). Assuming wi < W2, multiply the differential 
equation in (56 1 by € Hqq{I), then integrate by parts over /: 

1 



We deduce 11^+11/^(7) = 0, whence m < on / and finally Cwl{x) < Cw}{x) for all x e /, 
thus completing the proof of (i) . 

(3) Regarding (ii), it is customary to use the expression of c^'' given by Eq. (43 1. In particular, 
it is evident that we must have c^^(O) = 1 because Cw\o) = 1 for all w > due to the 
boundary conditions of Problem (42 1. Moreover, if a; G (0, 1] then 

c(„°)(a:) -e-'"^" {w ^ +oo), 

therefore Cw\x) — > for w +oo. 

(4) Fix now any cq G (0, 1). We observe that for w = we have c^\x) = 1 while, as a 
consequence of (ii), for w +oo it results Cw\^) 0. Since we know from (i) that 
w I— > Cw\l) is continuous, we conclude that Cw\i) takes all values in (0, 1] as w ranges in 
M+, hence there exists > such that Cw}{l) — cq. But (i) also tells us that w i-^ Cw\^) 
is nonincreasing, thus Cw\^) < cq for all w > w^,. Finally, since x Cw\x) is continuous 
on I with Cw\o) = 1 > cq, for all w > w^, there must exist a point x^, G / such that 
Cm "* (iu,) = Cq. Uniqueness of x^ follows from the strict monotonicity of Cw"^ (x) with respect 
to X G / (use Proposition [l] along with Cw\x) > for all a; G / to discover dxCw\x) < 
over I). This gives (iii). 

We observe that, under the hypothesis w > w^,, it results Cw\x) > Cq for x G [0, x^] 
and Cw\x) < cq for x G {xw, 1]. 

(5) Finally, let us consider statement (iv). First we check that the limit of x^, for w — > +oo 
exists by observing that for 'W2 > wi > it results 

•^Li ) — '^wl{^W2) ^ Cq ~ c[^^{Xwi) , 

where we have used the monotonicity of w i-^ Cw and the definition of 

tively. Since x i-^ Cw}{x) is nonincreasing, this says Xw2 ^ x^i, hence w i-^ is monotone 

nonincreasing and admits therefore a limit at +oo. Also notice that 

(57) lim Xw — inf a;^, G /. 

Assume now by contradiction that limi„^_(_oo Xj^j — > 0. Fix then w > w^^ and consider 

the point a; = ^ G /. Owing to (57l we have ^ < Xw, whence Cw\0 ^ c^^H^w) ~ co, for 
all w > , so 



c(^)(0= lim cW(O>co>0, 

w — *+oo 

which is in contrast with (ii). Thus — > and we obtain (iv). □ 



Theorem 9. Let Assumptions\^yjhold. Then there exists a positive solution to Eq. (46 1. 
Proof. Let us denote 

1 

Ciw)= fnc<^^)dx; 



since c'"q\x) = 1, we deduce C(0) = r(l) > in view of Assumption [i] 
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(1) We claim that C is continuous on M+. Indeed, let Wi, W2> 0. Using Assumption [s] and 
Proposition |8] we get 

1 

\C{W2)-C{w,)\< f |r(cL°))-r(cL"))|dx 



W2 ^Wi 

SO that for W2 Wi it results C{w2) C{wi) 



Due to the positivity of C(0), existence of a positive root of Eq. (46 1 is simply obtained by showing 
that C(w) < for w large enough: Continuity will do the rest. 

For Co e (0, 1) fixed by Assumption |4] let uj* > be the value of w mentioned by Proposition 
[Sj^iii): We prove that there exists w** > ui* such that C{w) < for all w > w**. To this end, we 
will henceforth assume w > Wt. and rewrite C as 

s„ 1 

CH= / Tic^^^)dx+ f T{c<^^)dx^:Ciiw) + C2iw). 



2™ 

(2) Notice that Ci(w) > for all w > w.,, indeed T{c^^\x)) > for a; S [0, x^). The Mean 
Value Theorem implies the existence of Xw G (0, x^) such that 

Ciiw) = x^T{cl°Hx^)). 

By consequence, using Proposition [sj^iv) we get Ci{'w) < TmXw for w +oo. 

(3) Concerning C2, we begin by observing that there exists w.^, > w■^, such that T{cw\x)) < 
for all X G [xw, 1] and all w > w*, with T{cw\x)) not identically zero on [xu,, !]• Indeed, 
we know from Proposition [sj^ii) that the c^^'s tend pointwise to zero almost everywhere 
in [0, 1] as w ^ +00, and from Assumption [3] that r(0) < 0. Consequently we also deduce 
C2{w) < for all w > iB^. 

(4) Next we claim that C2 is a nonincreasing function of w. Given W2 > wi > w^, we use the 
monotonicity of w 1— > Cw"^ and of F to discover indeed 

111 
C2iw,)= J nc(^l)dx> J F(cL"i)da;> J F(4°)) = C2(zz;2). 

X-(jj-^ iC-uj-j^ Xi(j2 

(5) Fix e > 0. Since Ci tends to zero for w +00, there exists w** > ui* such that Ci{w) < e 
for all w > w**. Moreover, C2(w**) < 0. Let us now choose e — ^|C2(w**)| = — ijC2(w**). 
If w > w,*, we find 

C{w) = Ciiw) + C2{w) < e + C2{w,,) = ]^C2{w,,) < 
and the thesis follows. □ 



7. Conclusions and perspectives 

In this paper we have developed a macroscopic model which describes the evolution of a tumor 
cord using the theory of mixtures. The tumor mass has been modeled as a three-phasic porous 
medium composed by cells, extracellular fluid, and extracellular matrix, growing along a blood 
vessel whence cells get the necessary nutrient for their vital functions. As a matter of fact, since 
the main interest was in investigating the growth process in connection with the availability of 
nutrient, the extracellular matrix has been regarded as a rigid non-remodeling scaffold in which 
cells live and proliferate and extracellular fluid flows, therefore it has not played a primary role 
in the dynamics of the system. By means of suitable principles of mixture theory, in particular 
the reference to a Darcy-like framework to model interactions among different components, the 
problem has been reduced to a system of two partial differential equations. The first, for the cell 
volume ratio, is derived from the cell mass balance in which the velocity is explicitly expressed 
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in terms of the internal stress state of the cells. It also includes a source/sink term describing 

proliferation or death of the cells according to the current size of the cell population itself, as 
well as to the available amount of oxygen present in the mixture. The second, accounting for the 
dynamics of the nutrient, is a linear parabolic equation modeling the diffusion of oxygen from 
the vessel wall toward the mixture and its contemporaneous uptake by the cells. This simple 
setting has to be understood in the light of the experimental evidence that the contribution of 
a possible transport by the extracellular fluid to the motility of oxygen molecules is definitely 
negligible with respect to their high diffusivity, hence the corresponding advection terms can be 
rightly dropped from the equations. Interaction of the growing tumor cord with a surrounding 
host tissue composed by normal non-proliferating cells has also been considered. 

The mathematical formalization of these ideas results in a free boundary problem constituted 
by a system of two partial differential equations supplemented by suitable boundary and interface 
conditions. The free boundary is represented in this context by the interface separating the cord 
from the host tissue, which changes in time according to the growth of the tumor in interaction 
with the external normal cells and has to be determined as a further unknown of the problem. 

The model is intrinsically multidimensional. In addition, it is not linked to a particular geometry 
because it does not assume any special configuration of the domain. In view of this, it is in principle 
applicable to a wide range of different scenarios. In this work, a two-dimensional problem has been 
specifically considered, focusing on the axial growth of the cord along the blood vessel and on its 
simultaneous expansion outward the vessel. Numerical simulations have shown ability of the 
model to reproduce specific features of vascular tumors, like the formation of an inner vital zone, 
constituted by sufficiently fed cells that duplicate, and an outer necrotic zone at the periphery 
of the cord, where conversely cells are not reached by a proper quantity of nutrient and cannot 
hence proliferate nor survive. It is worth pointing out that some different models of tumor growth 
assume a priori the existence of such zones and describe them invoking specific equations. In the 
present model, their formation is instead recovered a posteriori as a consequence of the overall 
dynamics predicted by a unique set of equations, by appealing to relatively general guidelines. 
Another interesting feature made evident by the numerical results concerns the evolution of the 
system toward a steady state. The cord clearly exhibits two different behaviors in the front part, 
the head, which is vital since always globally fed by a sufficient amount of oxygen, and in the 
rear part, the tail, which elongates as the head moves forward along the blood vessel, keeping a 
rectilinear shape with a nearly constant width in the transverse direction. 

These considerations suggest that, as far as axial growth is concerned, two main questions have 
to be addressed, namely the dynamics of the head in terms of shape and velocity of propagation and 
the equilibrium of the tail in terms of maximum penetration inside the host tissue. In the present 
work, we have mathematically investigated this second issue, studying existence and regularity 
of physically significant solutions that describe the distribution of cells and nutrient across the 
tail at the steady state. By addressing furthermore the free boundary problem via a perturbative 
expansion of the solution, we have also been able to characterize qualitatively and quantitatively 
the steady growth width of the cord. 

The model developed in this paper should be regarded as minimal, in the sense that it is 
deduced from few outline principles of a general well-coded theory, and takes into account only 
those really fundamental aspects to obtain a qualitative mathematical and biological description 
of the macroscopic phenomenon under consideration. Of course, it is liable to many improvements, 
which still at a basic level may involve, for instance, a more accurate description of both host cells 
and extracellular matrix, that here are essentially passive and lacking in their own physiology, 
and a more sophisticated coupling between the dynamics of the tumor cells and the nutrient, 
accounting for different uptake rates on the basis of the specific functions (survival, proliferation) 
carried out by cells. However, we believe that a minimal model, for which both a theoretical 
comprehension and validation are possible in view of its relatively simple structure, constitutes a 
solid foundation, hence a safe starting point, to tackle more complicated situations. 
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